Equation of State for a Complex Plasma from Recursive Bayesian Estimation 



(N 
O 

o 

Q: 

Q-r 



C/3 . 

in" 
'Hi 



o 

(N 



Neil P. OxtobyQ Elias J. Griffith, Celine Durniak, Jason F. Ralph, and Dmitry Samsonov 
Department of Electrical Engineering and Electronics, 
University of Liverpool, Liverpool, L69 3GJ, United Kingdom 
(Dated: 5 December 2012) 

We demonstrate tliat high-precision object tracfcing using recursive Bayesian estimation signifi- 
cantly improves empirical estimates for the equation of state of an experimental two-dimensional 
complex plasma, which we exposed to a series of electrostatic pulses to excite shock waves of different 
intensities. The estimate for the adiabatic index approaches that of a monatomic ideal gas (7 = 5/3), 
verifying a common assumption in the analysis of complex plasmas. Lower-precision dust kinematics 
from particle tracking velocimetry introduced a bias in the estimate of the adiabatic index. 
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An equation of state, such as the ideal gas law, is a 
mathematical relation between physical constants and 
macroscopically observable properties of a single phase 
of a system in equilibrium Equations of state are 
path-independent, and so can be explored by changing 
a system along any convenient intraphase path in state 
space between equilibria. Interphase paths include a 
phase transition — a discontinuous change in one or more 
system properties. For example, the significant increase 
in volume when liquid water evaporates. Nonequilibrium 
paths, whether intraphase or interphase, can also be used 
to infer an equation of state, but an assumption is re- 
quired to link the nonequilibrium states to the equilib- 
rium states. This is the case in shock wave physics where 
otherwise unreachable high pressure and high density re- 
gions of state space are explored. Pressure-density curves 
from shock experiments do not provide enough thermo- 
dynamic information to infer an equation of state Q , but 
can be used to fit parameters in an assumed equation of 
state. We explore parameter fitting in the ideal gas equa- 
tion of state, applied to a complex plasma. 

A laboratory complex plasma consists of plastic mi- 
crospheres suspended in a low-density ionized gas. The 
microspheres are often referred to as dust particles in 
analogy with dusty plasmas observed in astronomy 0, Q ■ 
Fast-moving electrons and relatively slow-moving ions in 
the plasma deposit a net negative charge on the dust, 
which repel each other via a screened Coulomb force 
(Yukawa or Debye-Hiickel) [Hj. Condensed- matter-like 
behaviour results when the dust is electrostatically con- 
fined, with the dust mimicking microscopic constituents 
of a fluid (individual molecules/atoms), yet being observ- 
able on a macroscopic scale (even to the naked eye) . The 
space between dust particles is occupied by a rarefied gas, 
so these dusty plasma structures experience very little 
damping, and are therefore considered to be representa- 
tive models of liquids and solids . Dusty plasmas are 
thus an excellent vehicle for exploring the microscopic 
kinematics of melting processes and crystal formation. 



Thermodynamic properties of the dust can be cal- 
culated from the kinematics of the individual parti- 
cles. Particle positions extracted from images determine 
the dust density via Voronoi analysis 0, Hj, andparti- 
cle velocities determine the kinetic temperature [3, ■ 
Dust kinematics are normally estimated using particle 
tracking velocimetry (PTV) where average velocity 
vvTw{t+T/2) = [x{t+T)-x{t)]/T is calculated from con- 
secutive position measurements x{t + T), x{t), which are 
extracted from a sequence of images taken with a high- 
speed camera at a frame rate of l/T (typically 500-1000 
frames per second). The velocity calculated in this way 
is subject to two sources of inaccuracy: position uncer- 
tainty in the measurement, and nonzero acceleration. For 
very high frame rates T — 0, upTV is limited by position 
uncertainty, which is due to finite pixel size and noise in 
the camera sensor [l2, 13|. These limitations can lead to 
artefacts in results calculated from PTV-estimated kine- 
matics. Recursive state estimation (also known as ob- 
ject tracking Hlli) has been employed to estimate the 
kinematics of complex plasma particles [i, [3. Object 
tracking algorithms filter noisy measurements via a set 
of equations to produce estimates of the instantaneous 
kinematics which are resilient to the limitations discussed 
above. The most ubiquitous recursive Bayesian estimator 
is the Kalman filter [l3|. 

In this work we employ object tracking using an in- 
teracting multiple model tracker based on Kalman 
filtering (KF) to generate thousands of simultaneous 
particle tracks from shock wave experiments on a two- 
dimensional (2D) complex plasma. We use Rankine- 
Hugoniot relations to calculate Hugoniot curves aris- 
ing from the estimated kinematics, and interpret the re- 
sults using the ideal gas equation of state. We compare 
our KF results with those from PTV, which does not ac- 
count for the non-negligible particle acceleration in shock 
wave experiments, and find that PTV introduces a sys- 
tematic error that gives rise to a bias in the estimated 
equation of state. Object tracking can incorporate parti- 
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cle acceleration in one of two ways: either by modelling 
known interparticle forces, or by including acceleration 
in the recursive estimation algorithm along with position 
and velocity. Multiple forces govern the dust dynam- 
ics , including Yukawa (particle interactions) , external 
(confinement and excitation) , and friction (plasma drag) 
forces. For a well-characterized system with known par- 
ticle charge, Debye length, damping rate, particle mass, 
etc., the total force on each particle is calculable from 
the positions and velocities of all particles. Such known 
accelerations can be exploited by a state estimation fil- 
ter to improve the accuracy of particle tracks 0. How- 
ever, complete characterization of a complex plasma is 
extremely challenging. Fortunately, even in the absence 
of exact knowledge of the physical parameters, object 
tracking allows the acceleration of the particles to be es- 
timated recursively — along with position and velocity 
— thereby improving force estimates. This is the situa- 
tion considered here. 

Equation of State — The most common thermody- 
namic relation for compressible flow is the ideal gas law, 
which can be written in terms of the following specific 
(per unit mass) quantities: pressure p, internal energy e 
and density n as 



(1) 



where 7 is the adiabatic index. In the case of an ideal 
gas undergoing an isentropic process (reversible adia- 
batic process), 7 is the ratio of specific heats at con- 
stant pressure and constant volume. Deviations from 
the ideal gas law were first considered by van der Waals 
in 1873 to account for finite particle size and molecular 
interactions . However, the ideal gas law can be ap- 
plied to non-ideal gases with sufficient accuracy in many 
cases [iOl- 

Here we explore the p(e,n) relation in a non- 
perturbative manner by exciting a series of shocks of dif- 
ferent magnitudes in the dust O HH, A normal 
shock wave is one where the shock front is normal to 
the direction of propagation, and the bulk flow is one- 
dimensional. In the frame of a normal shock moving 
at speed us, the Rankinc-Hugoniot shock relations for 
conservation of mass, momentum and energy across the 
shock front arc, respectively, [31 
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where ui/2 ~ ug — ui/2 denotes the particle speed down- 
stream/upstream (ahead/behind) from the shock in the 
laboratory frame and the moving frame is denoted with a 
tilde as Ui/2- Speciflc density n and internal energy e are 
equal in the laboratory frame and the moving frame, but 
pressure has a kinetic component. We found this to be 



small (<J>%) compared to the Yukawa pressure (shear 
stress [23[) so that p p. Equation ((2c)) is known as 
the Hugoniot [23 |. Using equation ([1]) to eliminate in- 
ternal energy from equation (|2cl). and combining with 
equation (j2bp , we can write UJ, 



V (7 + 1) - (7 - 1) 

(7 + 1) - 77(7 - 1) 



(3) 



where ^ = P2/P1 is the shock strength and rj = n2/ni is 
the compression ratio across the shock. A least-squares 
flt of the experimental results to equation ^ yields an 
estimate of 7, and hence an equation of state for the 
shocked dust in the form of equation ([1]). 

Polytropic index — With the ideal gas law as the 
equation of state, the polytropic index g can be used to 
describe the physical nature of a process that changes an 
initial state (downstream pi,ni) to the flnal state (up- 
stream P2,n2). Polytropic processes follow p/n^ = C [JIj 
which describes a curve in the pressure-density diagram, 
with g and constant C defining a solution for the pres- 
sure/density changes linking the initial and final states. 
Equating initial and final states (both equal to C) then 
combining with ^ and 77 then solving for g allows the 
polytropic index to be expressed in terms of ^ and 77 as 
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where g = indicates an isobaric process, 17 = 1 is an 
isothermal process and = 7 is an adiabatic process. 

Experiment — The experiment involves levitating a 
2D array of microspheres 10mm above the floor of an 
Argon-filled chamber pressurized to 2.05Pa. The spheres 
(9.2/xm diameter) were then allowed to settle into a well- 
spaced crystalline structure, forming a "plasma crys- 
tal" (2^ which is visible to the naked eye when illumi- 
nated by a laser sheet (figure [T]) . These dust particles 
each hold an approximate charge ofQ = IGOOOe and have 
a Debye length of Ajj = 1.0mm [2^ [2^. Shock waves 
were created by an electrode located to one side of the 
visible area, which was pulsed for 2 seconds with a volt- 
age selected from —20V to —50V in 5V steps. The crystal 
was allowed to reset between each run (requiring approx- 
imately 100 seconds). The experiment was repeated at 
each voltage level to reduce the impact of local variation 
in crystal structure that can form on reset. The dust was 
imaged from above by a grayscale camera at 500 frames 
per second for 1.2 seconds, and the resulting images pro- 
cessed by PTV and our Kalman-filter-based tracking al- 
gorithm to obtain the dust kinematics. Further details 
of the experimental apparatus required to achieve this is 
described in references 



22| and [27[. An example snap- 



shot image of the dust, enhanced for presentation with 
enlarged dots and false color, is shown in figure [T] with 
a zoomed inset around the shock front. The symmetry 
inherent in normal shocks permits a ID description of 




FIG. 1. Typical overhead image witii zoomed inset, enhanced 
for presentation (enlarged dots, false color). The visible area 
is 32.8mm square imaged at a resolution of 1024x1024 pixels. 



the thermodynamics. The profile values were calculated 
as robust average quantities (median, median absolute 
deviation) in each of 50 bins equally spaced along the X 
axis, which spanned the Y axis. The density and pressure 
profiles were used to generate a p{n) plot below. Density 
is the inverse of the Voronoi cell area 0, Q , and pressure 
is calculated as shear stress, i.e., the X-component of the 
2D stress tensor (23{ . 

Our investigation proceeded as follows. The shock 
front was identified as a peak in the density profile evo- 
lution (see figure [2]), from which the shock position and 
speed was determined. The upstream and downstream 
quantities in the Rankine-Hugoniot shock-jump relations 
were selected from 0.656mm (1 bin) behind the shock 
front and 3.28mm (5 bins) ahead. We needed to look 
further ahead to overcome the finite width of the shock 
front (an ideal shock would have vanishing width). Re- 
sults were also sensitive to the chosen upstream distance 
due to structure a few millimetres behind the shock (see 
multi-shock discussion below). 

A principal Hugoniot p{n) requires repeated shock ex- 
periments of different magnitudes, sharing a common ini- 
tial condition. Reliably reproducing the same initial con- 
dition for a dust crystal is extremely difficult, if not im- 
possible. For this reason the data was post-selected from 
the densest cluster of initial conditions and constrained to 
lie within 1% of the cluster centroid. This is illustrated in 
figure |3] where the post-selected downstream data points 
are shown as blue dots and all others as red crosses. The 
inset of figure [3] also shows the post-selected upstream 
data. From 13 similar experimental runs, 118 total data 
points were generated, of which 26 were post-selected. 

Results — Figure [5] clearly shows two density peaks: 
one leading (black) and one trailing (white). Such multi- 



FIG. 2. Evolution of the dust particle number density profile 
n{X,t). The density peaks (crosses) clearly shows the prop- 
agating shock wave (black), with an additional trailing shock 
(white). Quadratic least-squares fits are shown. 
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FIG. 3. Initial pressure and density (ni,pi) for each run: 
blue dots survived post-selection (see text). Inset: Pressure- 
density diagram showing all post-selected results: initial (blue 
dots) and final states (black circles). 



shock structures [28| can be described by a sequence of 
jump relations of the same form as ([2]), and are stable 
when each wave travels at least as fast as the trailing 
waves (so the lead shock is not overtaken). This was the 
case here, where least-squares fits for the shock (S) and 
trailing (T) peak positions determined the wave speeds 
(mm/s) to be us(0 = -17.2t-|-43.7 and UT{t) = -8At + 
33.0. The shock front is detectable after 0.24 seconds, 
with the trailing shock resolvable after 0.45 seconds. 

Figure m shows shock strength vs. compression for both 
PTV and KF (object tracking). Least-squares fits to 
equation ([3]) determined the estimated adiabatic index 
for each approach. We found 7kf = 1-67 ± 0.01, 
which is consistent with that of a monatomic ideal gas 
7idcai = 5/3 = 1.66. By contrast, PTV overestimated 
the value at 7ptv = 1-79 ± 0.01. The following discus- 
sion explores the reason why PTV overestimates 7 from 
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the ^-r] plot. The post-selection renders the downstream 
data not culpable since variation is minimal by design. 
So, ^ and rj are directly proportional to P2 and 712, respec- 
tively. Dust pressure (and hence £,) is dominated by the 
Yukawa interaction, which is non-linear in interparticle 
spacing r (e.g., see [§]), and particularly sensitive to er- 
rors in r. When averaged, these errors propagate through 
the non-linearities to create a biased overestimate for the 
pressure, thereby shifting erroneous results upward in the 
^-77 plane. Dust density (and hence 77) is easily underes- 
timated when excited dust particles intermittently leave 
the plane of illumination. The object tracking process 
provides a robust way to maintain tracks for these parti- 
cles that PTV does not. The reduction of 77 due to such 
track loss shifts erroneous results to the left in the ^-ry 
plane. In figure S] the PTV result is above and left of the 
KF result (which may lie above or on the true result). 

We determined the polytropic index of the shocked 
dust from our KF results via the mean value of equa- 
tion (g]). We found ^kf = 1-71 ± 0.07 which satisfies 
(7 « 7, thereby demonstrating that shocks in a dusty 
plasma (which here behaves as an ideal gas) constitute 
an adiabatic process. Thus, our investigation experimen- 
tally verifies the expectation that shocks in an ideal gas 
are an adiabatic process Conventionally this is often 
assumed to hold, at least as an approximation, whereas 
we have verified the validity of this assumption experi- 
mentally. 
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FIG. 4. Shock strength vs. compression ratio for PTV and 
KF. Least-squares fits to equation ([3} give the adiabatic index 
7, with 3(7 confidence regions in each case shown in gray. 



Figure [5] shows the shock Hugoniot 2^ S^] , where the 
shock velocity us is linearly related to the upstream ve- 
locity U2 behind the shock front: Ms = Su2 + Cq. Here 
Co is the zero-pressure bulk speed of sound (speed of 
sound in an unshocked sample), and S* is a dimensionless 
constant of proportionality for the linear fit. The PTV 
and KF methods estimate Cq to be 26.8 mm/s and 21.8 
mm/s, respectively. These values are in line with the 25 
mm/s and 28 mm/s speeds of sound observed in [sij and 



[32l | via different techniques. The very low speeds are 
due to the extreme softness of the dust crystal, which is 
due to the large interparticle spacing in comparison with 
the particle size [sij. (Recall that the particles shown 
in figured] have been enlarged for clarity.) In the fits of 
figure [SJ the coefficient of determination showed that 
the KF data (i?^ = 0.64) followed the expected linear 
trend far better than the relatively scattered PTV data 
(i?2 = 0.27). While there is stih some way to go (i?^ = 1 
is a perfect fit), this result further increases our confi- 
dence in using object tracking methods to analyse dusty 
plasma experiments, rather than the prevailing PTV ap- 
proach. 
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FIG. 5. Shock front speed vs. upstream particle speed {Shock 
Hugomot) for both PTV and KF (object tracking). 

Conclusion — We performed shock experiments in 
a two-dimensional dusty plasma. This model system 
for exploring the microscopic dynamics of molecular sys- 
tems allowed us to calculate thermodynamic variables for 
the dust (pressure, density) directly from the kinemat- 
ics of the myriad dust particles. The kinematics were 
estimated using two techniques: a recursive Bayesian 
state estimation technique (object tracking), and particle 
tracking velocimetry (the standard approach in complex 
plasma physics). Conservation laws known as Rankine- 
Hugoniot equations were combined with the ideal gas law 
to estimate the adiabatic index of the dust, which veri- 
fied two common assumptions: that a dusty plasma can 
be considered as an ideal gas, and that a shock wave in 
an ideal gas is an adiabatic process. 

We found that PTV analysis overestimated the adia- 
batic index (and hence produced an inaccurate estimate 
for the equation of state), due to a systematic error intro- 
duced when the errors from PTV are used in the nonlin- 
ear equations. This reinforces that object tracking pro- 
duces reliable dust particle tracks @ and should be used 
when analyzing complex plasma physics experiments. 
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